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Computation of nucleation of a non-equilbrum first-order pliase transition using a 

rare-event algorithm 

David A. Adams'^ Robert M. Ziff^'Q and Leonard M. Sander°[t] 
"^Department of Physics and ^Department of Chemical Engineering, 
University of Michigan, Ann Arbor MI 48109-2136. 

We introduce a new Forward-Flux Sampling in Time (FFST) algorithm to efficiently measure tran- 
sition times in rare-event processes in non-equilibrium systems, and apply it to study the first-order 
(discontinuous) kinetic transition in the Ziff-Gulari-Barshad model of catalytic surface reaction. The 
average time for the transition to take place, as well as both the spinodal and transition points, are 
clearly found by this method. 



I. INTRODUCTION 



In many systems, rare events occur with a very low probability compared to typical events. Sometimes they are 
of central interest. Examples include the extinction of diseases [1] or of populations [2], network queue overflow 
I— —I [3], and slow chemical reactions [1]. The study of such processes poses a particular challenge to simulation. In 
i-G the field of chemical physics, many rare-event techniques are commonly used: transition path sampling 5 , transition 

^ interface sampling [6j , milestoning [7] , the string method |8] , and the weighted-ensemble method [9^ , to highlight a few 
methods. A thorough review can be found in [10]. Most of these methods require that the system being studied has 
an underlying energy landscape, which precludes their use on non-equilibrium systems, i.e., systems that lack detailed 
balance. Forward Flux Sampling [11 (FFS) is a rare-event technique designed specifically for non-equilibrium systems, 
and has proven useful in studying genetic switches [m - fT3] , nucleation p"3l - fT5] , isomerization of alanine dipeptide [16] , 
^ and the Maier-Stein model of reaction dynamics [13l|T7l[l8]. In this paper, we develop a variant on FFS and use it 
to study the first-order non-equlibrium phase transition in a catalysis model. 

FFS was developed to measure transition rates between two locally stable regions A and B separated by a high, 

C featureless barrier. When the barrier separating these regions contains long-lived metastable states, FFS is generally 

r^ inaccurate or inefficient, depending on how it is applied. We have overcome this limitation of FFS with a variant, 

(~| which we call Forward Flux Sampling in Time (FFST). In our method, we adjust for long-lived metastable states by 

O measuring the times associated with sampling the region between A and B in the second stage of the FFS algorithm. 

^ O^ The method is described in detail in Appendix A. 

We apply FFST to the Ziff-Gulari-Barshad (ZGB) catalytic surface-reaction model [T^. This model is of interest 

T-H because it has a first-order phase transition which acts in many ways like an equilibrium phase transition (it shows 

^ critical behavior, nucleation, etc.), but the model is manifestly non-equilibrium. Many techniques have been applied 

^"^ to the transition in order to tease out its properties. FFST allows us to study the dynamics instead of overall rates. 

^^ Using FFST we found transition times for nucleation as large as 10''° Monte Carlo steps (MCS), more than 30 orders 

of magnitude longer than those accessible to direct simulation. The method generates not only the transition times 

but the ensemble of most-likely states as the system progress from one phase to another. This allows us to measure 

^^ properties of the ensemble during the transition, which helps determine the pathway. 

The outline of this paper is as follows: In section |Tl] we describe the model and simulation method. In section \Ln\ 
we describe our results, and in section |IV| we summarize our conclusions. In the Appendix A we give the details of 
the FFST technique, and in Appendix B we test FFST on an exactly soluble one-dimensional system. 
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II. MODEL AND SIMULATION METHOD 



The ZGB Model 



The ZGB model was introduced to study the behavior of the oxidation of carbon monoxide (CO) on platinum 
surfaces (TH]. The catalytic surface is represented by a square lattice on which CO and O2 can adsorb. The CO takes 
up one lattice site whereas the O2 dissociates into two O atoms which take two adjacent vacant sites. Both species, 
CO and O, are bound to the surface until the other species adsorbs on a neighboring site. At this point the CO and 
O form CO2 and desorb from the catalyst leaving two lattice sites empty. The rates of reaction CO + O ^ CO2 and 
the desorption of CO2 are assumed to be infinite. The state of the catalyst is controlled by the fraction of the time 
CO is attempted to be placed on the catalyst; this fraction is called pco- 

In this basic model, there exists a region of steady-state reaction, bordered from below at pco = Pi by a second- 
order, continuous kinetic phase transition to an 0-covered state, and above pco = P2 by a discontinuous, first-order 
kinetic phase transition to a CO-covered state. The first-order transition is robust to small changes in the model 
(diffusion, small desorption, etc.) and is seen experimentally at low temperatures as a sharp transition from high 
to low reactivity [201. The second-order 0-poisoning transition is weak and not seen experimentally. The first-order 
transition is associated with many complex oscillatory and wave phenomena in theoretical [21, 22; and experimental 
systems [531 HI] and has thus received much attention. It also serves as a paradigm for general first-order kinetic 
phase transitions |25j . 

Studying the first-order transition has proved to be a challenging problem in simulations. Because of the difficulty 
of nucleating a sufficiently large CO cluster or island, simply increasing pco from the reactive steady-state misses 
the transition point, and instead the CO-poisoning is seen to occur at pco ~ 0.5277 [2S]. That point is close to an 
effective spinodal point p* , where the transition occurs without any kind of barrier [27j [28] . 
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FIG. 1: The poisoning time T vs. pco for various values 
of system size L and FFST trials A''. 



FIG. 2: The poisoning time T vs. pco for various values 
oi LwithN = 10^ 



Estimating the value ofp2 and properties of the first-order phase transition have received significant attention. Using 
a "constant-coverage" technique, the values p2 = 0.52560(1) [571 HI] and more recently 0.525615(5) [2S] have been 
found. Other methods, including histograms |30j . epidemic analysis [31j . and epidemics and trigger waves j32l I33j . 
have also been used to probe the first-order transition. 

The results of these studies is that there is a first-order transition at p2 ~ 0.5256, and a spinodal at p* ~ 0.527—0.528 
[IHKin], although the precise, and somewhat lower, value p* = 0.52675(5) has also been proposed [31]. For p2 < pco < 
p* , there presumably exists a critical CO-cluster size, below which clusters tend to shrink and above which they tend 
to grow. That critical size changes from 00 to something of the order of the lattice spacing as pqq goes from p2 to 
p* . For finite systems, the behavior is controlled strongly by the boundaries. Using periodic boundary conditions, 
as Pco is increased, the largest CO-cluster goes from being isolated, to wrapping around one direction (leading to 
two interfaces that are fiat on the average) to wrapping around in both directions, as illustrated in Fig. ^ In the 
intermediate coverage region, the constant-coverage method gives p2 accurately with very small finite-size effects. 
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FIG. 3: The scaled transition time InT/L vs. L for various values of pco- 



While the constant-coverage technique maps out the transition, it does not provide any information about the 
dynamics of the system. For that, it is necessary to study the standard (constant-rate) ensemble. But in that case, 
the nucleation barrier makes it virtually impossible to study dynamics except for very close to the spinodal point. To 
overcome this challenge, we use a modified FFS technique to find nucleation dynamics as well as overall rates for a 
wide range of pco values. 

B. The simulation method 

We study the ZGB model using the constant-rate ensemble on a square L x L lattice with periodic boundary 
conditions. The dynamics involve repeated attempts to adsorb the CO or O2 species. The procedure for an adsorption 
trial is given below. 

• pick r G [0, 1), if r < pcOi attempt to place a CO molecule, otherwise attempt to place O2. 

• pick a random lattice site {x,y) and continue if the site is empty. If placing O2, also pick a neighboring site 
(x ± 1, 0) or (x, y ± 1) and continue if that site is also empty. 

• Place the CO or O2 (dissociated) onto the empty lattice site(s). 

• For each lattice site now occupied, check all of the neighbors of that site and determine if any of those neighbors 
are of the opposite species. If any are, remove a randomly chosen neighbor of the opposite species with the 
recently placed species. 

Before any adsorption trials can be attempted, the system must be initialized. Because we wish to study the 
first-order phase transition, we prepare the system in the reactive state. Starting from an initially empty lattice 
with pco close to the spinodal would frequently lead directly to the poisoned state. Instead, we prepare the system 
by adsorbing CO with probability Pqq = 0.07 and O with probability Pq = 0.43 on every site. This initialization 
generates invalid states with CO neighboring O. To remedy this problem, we run the simulation for 10 MCS, which 
drives the system to a valid state. These 10 MCS "burn off' a significant fraction of the original CO and O including 
those with neighbors of the opposite species. 



C. Forward Flux Sampling in Time 



To study the first-order phase transition, we use the FFST algorithm, which is described in detail in Appendix A. A 
sketch of the algorithm is as follows. Before the simulations begin, we define a starting region in state-space A and an 
ending region B, bounded by 'barriers' Aq and Am, respectively. We also define Ai through Am-i as dividing surfaces 
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FIG. 4: The average largest cluster size M vs. coverage C 
for various values of pco and L — 96. 



FIG. 5: The spanning probability vs. C for various values 
of pco and L — 96. 



in state space that effectively mark the distance between A and B. The first step of the algorithm is to run a long 
simulation starting in A and recording where the sample path crosses Aq going out of A, and the average time spent 
after crossing back into A before leaving again; we call this the internal return time Tint- Next, we start sample paths 
along Aq where the initial simulation crossed, and continue them until they reach Ai or go back into A. The fraction 
of paths that reach Ai gives an estimate of the probability of reaching Ai without going back into A, P(Ai|Ao). We 
also keep track on the average time it took to reach Ai and go back into A starting from Ai. In the second step, we 
continue paths from the locations along Ai, where the previous paths stopped, and run them until they reach A2 or 
go back inside A, crossing Aq. The results give estimates for P(A2|Ai) and the time it takes to reach A2 or A from Ai. 
This process is repeated, step by step, until Xm is reached on the M"^ step. Finally, we use the results collected to 
calculate the overall transition time, which is given by the probability of reaching Am from Aq without going into A 
(P(Am|Ao)) times the average time it takes to leave A. The overall transition probability is given by the product of 
the intermediate transition probabilities, P(Am|Ao) = Ili^o P{^i+i\^i)- The time to leave A is the average time it 
takes to return to Aq from inside A, T^^t^ plus the average time it takes to return to A from outside Aq, Tf,,j.t, which 
we calculate from the times measured during the second stage of the algorithm. The gains of FFST over FFS are 
illustrated in Appendix B. 

Thus, the technique follows fruitful paths from the reactive state to the CO-poisoned state. To apply the FFST to 
the ZGB model, we must first define an order parameter, which is used to determine progress towards the poisoned 
state. This function should smoothly increase as system transitions from being reactive to poisoned. We chose the 
fraction of CO on the lattice (C) as our order parameter. 

FFST uses the order parameter to make surfaces (barriers) in state-space which are used to mark progress. Simu- 
lations are run from locations along a barrier in state space until they reach the next barrier, as defined by the order 
parameter, or the first barrier. 

These barriers can be placed before the sampling begins, which we call static barriers. The disadvantage of static 
barriers is that without a priori information about how to best place them, some barriers will have a large effective 
separation creating a performance bottleneck. To overcome these performance issues, two methods have recently 
been developed [351 136| which determine where to place the barriers automatically. To measure the poisoning time, 
we use the dynamic barrier placement method from |36] to place the barriers during the FFST algorithm. We use 
static barriers for measurements of the average largest cluster, largest cluster spanning probability, and committor 
probability fS] ; these terms are defined below. We use static barriers for the obserables because we want measurements 
at uniformly spaced intervals of coverage. When placing the barriers dynamically, we space them such that typically 
10% of trials make it to the next barrier before returning to the steady-state of the reactive region with low CO 
coverage. We use N trials per step in the FFST algorithm, except for the first barrier, which we use lOA^; we use 
7V/10 trials to determine where the next barrier should be placed. We locate the first barrier at the largest value of 
the order parameter found in the first 50 MCS. For the calculation of poisoning times, we used N ~ 10'^ and 10''. 
N — 10** was used for average largest cluster, largest cluster spanning probability, and committor probability. 
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FIG. 6: The average largest cluster size M^ vs. coverage 
C, as a function of L for pco = 0.5256. 



FIG. 7: The spanning probability ps vs. coverage C for 
various L and pco = 0.5256. 



III. RESULTS 



We desire the transition times from the reactive to poisoned state, which we call the poisoning time T. We determine 
T for a range of L: 32, 48, 64, and 96, as well a number of values of pco in the critical region, p2 < pco < P* , and in 
the reactive region, pco < P2 down to pco = 0.505. The results are shown in Figs. [T] and [2] 

Figs. [T] and [2] show that the poisoning time grows smoothly as pco decreases, with no indication of a transition 
at p2- Even ii pco is much smaller than p2, we found it possible to measure the poisoning time, demonstrating that 
the reactive state for finite systems is always metastable. We were able to determine maximum poisoning times of 
T — 10'^° - lO""^ MCS. The poisoning times converge to a value that is independent of system size for pco ~ 0.5275, as 
shown in Fig. [2] We associate this point with the spinodal point p* : as it represents the spontaneous and simultaneous 
nucleation of multiple clusters throughout the system, and is not influenced by the boundaries of the system. 

We attempt to find the form of the transition time as a function of L and pco- In niany non-equilibrium systems, 
transition times take the form T sa e^^ , where W is an effective energy barrier and L characterizes the size of the 
system. Fig. [s] shows InT/L which is an estimate for W. We see that InT/L does not have a strong system size 
dependence near p2, which may indicate that there is an effective Arrhenius energy. If this were the case, then InT/i 
will become independent of L for very large system sizes. 

Although FFST is primarily a tool for obtaining transition times, it also gives ensembles of states along each 
barrier as it progresses towards poisoning. These ensembles at a fixed value of coverage represent essentially what 
the constant-coverage ensemble attempts to mimic. Expectation values of quantities, like largest cluster size, can be 
taken over these ensembles, which give insight into the dynamics of the phase transition. 

We ran FFST with evenly spaced static barriers and measured the average largest cluster size and the spanning 
probability for the ensemble captured on every barrier. Spanning occurs when the largest cluster wraps around the 
periodic boundary and touches itself. We see a dependence on pco and L upon the spanning probability, as shown 
in Figs. [5] and [7] Fig. |5] shows that the smaller the value of pcO: the more likely that the largest cluster will wrap 
earlier in the path to poisoning, at fixed L. Fig. [7] shows that increasing L appears to narrow the range of coverage for 
which spanning has a non- negligible probability of occurring or not occurring, at fixed pco ~ P2- This suggests that 
the variation in shape of the largest cluster decreases with L. It also appears that reaching a spanning probability of 
50% is achieved at C « 0.42 independent of L at p2. This indicates that the average largest cluster is significantly 
elongated when it begins to span the system; as seen in Fig. |8] We also found a linear relationship between the scaled 
average largest cluster size and coverage, independent of both L and pco, as shown in Figs. [2] and |6] 

The ensembles of states at different values of the coverage can also be used to directly measure progress towards 
poisoning. By running every state in the ensemble until it returns to the reactive {A) or poisoned (B) state, one 
obtains the probability of poisoning from these particular values of the coverage. This probability as a function of 
the order parameter is called the committor probability, pb- (In the case of the ZGB model, this is only an estimate 
because the model is non-equilibrium, so the forward-tending ensemble obtained is not necessarily the same as the 
steady-state non-equilibrium distribution along the barriers.) We measured the committor probability for various L 
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FIG. 8: The evolution of the most hkely path to the poisoned state for L = 128 and pco = 0.5256 « p2- The pictures are in 
increasing values of coverage fraction from left to right starting at 0.09 and increasing in steps of 0.03. 

and pcOj a-s shown in Fig. l9] We find that for pco below p2, the larger the system, the larger the coverage must be 
in order to have a particular probability of poisoning. The opposite effect is found above p2 ■ At p2, the committor 
probability tends to a single form for the largest system sizes. This shows the significance of P2i as the point where 
Pb is independent of L. 

By keeping track of which state on a given barrier is responsible for a particular state on the next barrier, we were 
able to piece together complete trajectories from the reactive state to the poisoned state. Among these reconstructed 
paths it is possible to determine a most likely path. Every state ji on every barrier A^ has some probability of 
continuing to the next barrier before returning to the reactive state denoted by P-'(Ai+i|Ai), which is measured 
during FFST. We define the most likely path as the path which connects Aq and Am, which has the largest value of 
the product of the intermediate barrier crossing probabilities. Hi ^"'('^i+il'^i)- The states in Fig. ^show the most 
likely path for L — 128 at the transition point p2- The preferred path to poisoning involves wrapping around the 
system, then expanding to complete CO coverage. This path is favored because nucleating droplets have an effective 
kinetic surface tension which causes droplets to be unfavored. Once a droplet spans the system and the net curvature 
disappears, the cluster is significantly more favored. This preferred pathway of wrapping and expanding has also been 
seen in magnetic memory switching ^7\. This behavior is to be expected, at least near C sa 0.5, where the cluster 
usually wraps around in one direction. In that case, at the transition point p2 , the system should be equally likely to 
poison or return to the reactive state, so pb ~ 0.5 independent of L. 



IV. CONCLUSIONS 



In this article, we introduce forward fiux sampling in time (FFST). We use it to analyze the first-order phase 
transition in the ZGB model. We found a size-independent poisoning time at p* w 0.5275, which is associated with 
the spinodal point. The poisoning time is a continuous function of pco near the first-order transition. By inspecting 
the ensembles of states measured at each barrier, we found a linear relationship between the scaled average largest 
cluster size and the total CO coverage, which is practically independent of L and Pco- When the paths from barrier 
to barrier are connected, they make the ensemble of successful trajectories. We found the most probable of these 
trajectories and found that wrapping and then expanding is the preferred path to poisoning. 

At Pco — 0.5256, the committor probabilities appears to be independent of system size for large L. We believe 
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FIG. 9: The committor probability ps vs. C for four values of pco and different L. N = W^ for all data shown. 

that this is a signature of a first-order transition for the foUowing reasons; first, the committor probability is 0.5 for 
half coverage, which is has been previously used to determine the transition point |19j . Second, the matching of the 
entire curves for various L stems from the critical cluster size being infinite at the transition point. This implies that 
the probability of a droplet growing is always less than 1/2. The growth probability cannot be strongly dependent 
on the cluster size as clusters of arbitrarily large size must all have roughly the same growth probability, slightly less 
than 1/2. For large lattice sizes, what matters is how close the cluster is to spanning, which is only dependent on the 
mass of the largest cluster scaled by the total system size. Thus, we have found that the committor probability can 
be used to locate the first-order transition point of the system. Lastly, we found size-dependent poisoning times for 
systems well below the transition point. 

For the ZGB model, we found that the efficient implementation of FFS, i.e., a fixed number of crossings in the 
first stage, gives accurate results. This shows that ZGB doesn't have extremely long-lived metastable states, which 
is the primary advantage of FFST. But, we still found that FFST can outperform FFS in terms of smaller variance, 
which translates into better computational efficiency to reach a target variance. Specifically, we found that FFST was 
effectively 35% more efficient than FFS for a test case. 
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VI. APPENDIX A 



THE FFST ALGORITHM 



The FFS algorithm works well for problems with a featureless barrier, i.e., problems with infrequent but fast 
transitions. But, in the case of problems with slow transitions, typically caused by long-lived metastable states, it 
has significant shortcomings [3S1. In these cases, FFS either grossly underestimates the transition time or becomes 
nearly as inefficient as direct simulation (see Appendix B for an example). In this case FFS produces a constant gain 
in simulation efficiency, of the order of ten times faster, whereas the 'efficient' approach, where a fixed number of 
crossings is used to calculate the flux rate, can be 10^*^ times faster, as we found in this paper. In this Appendix, we 
derive the version of FFS, which we call FFST. We then outline two other variants on FFS which helps us illustrate 
the connection between FFS and the barrier method which we recently introduced [36] . 

As we explained above, FFS separates the problem of finding rare transition events into two problems: First, 
finding the rate of leaving the initial region A, and second, finding the probability of reaching the final region B from 
the surface of A without going back into A. Finding the rate of leaving A typically involves running single a long 
simulation until it has exited A a fixed number of times, say lOA'^ where N is the number of trials per barrier in 
the second step. The number of crossings divided by the total time of the simulation (discounting paths that reach 
B) gives an estimate of the rate of leaving A. This calculation of the flux is accurate only if the trajectories sample 
capture the important times and features of the entire landscape. For example, if none of the trajectories generated 
sample a long lived, metastable state between A and B, the estimated flux would be higher than the true flux. We 
show such an example in Appendix B. 

A simple and correct way to avoid this problem is to run the initial simulation until it reaches B one or more 
times. However, this is impractical because it effectively solves the problem using brute force, which is what FFS was 
designed to avoid, and would lead to miniscule efficiency gains over brute-force simulations. Our version of FFS has 
the accuracy of the 'correct' calculation and the efficiency of the usual ffux calculation. 



A. Forward Flux Sampling in Time 



Forward ffux sampling in time (FFST), the algorithm we use in this paper, performs FFS in terms of transition 
times instead of rates. The problem is decomposed in the same way as FFS[TT], utilizing the idea of endpoint regions 
A and B. In this view, there are three important times: the time to return inside of A from the surface of A (Text), the 
time to get back to surface of A from just inside A {Tint), and the time to reach B from the surface of A (Tf) without 
going back into A. The only additional quantity needed to calculate the transition time, Ttot, is the probability p of 
reaching B from the surface of A. Then: 



Ttot — ( 1 ) (Text + Tint) + Tf. 



(1) 



9 

If the first term, which represents the total time spent on unsuccessful attempts to reach B, is much larger than the time 
of a successful attempt to reach B, which is explicitly used in the construction of FFS, [(1/p) — l]{Text + Tint) ^ Tf, 
and p, is small then the above formula simplifies to the inverse of the FFS formula, 

Ttot ~ -{Text + Tint) = TJTT p^-^ ^. = 7 ■ (2) 

p P(Am|Ao) y^AfiJ KAB 

The efficiency gain of FFST comes from measuring Tint during the first step, and measuring T^xt a-nd p in the 
second step. Measuring Tint involves counting the time that the initial simulation spent in region A. In the second 
step, paths are allowed to run from the i"^ barrier to the (i + 1)**^ or back to A. We record the times to make these 
transitions to the next barrier (t*^^) and back into A (t°) in addition to the probability of making it to the next 
barrier P{\i^i\Xi). Text can be calculated by properly weighting the time it takes for different paths to return to A, 
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(3) 



1- n P(A,+i|A.) 



The numerator in the last expression is the sum over return times of all paths sorted by largest excursion. The 
probability that a path will start at Ao, make it to A^ and then return to Aq without making any further progress 
is product of the probability of reaching A^ {Y[j=q Pi^j+i\^j)) ^^'^ the probability of then returning to Aq without 

reaching A^+i (1 — P(Ai+i|Ai)). The average time of this excursion is the time it takes to reach A.; (X]/c=o ^k^^) P^^^ ^^^ 
time it takes to return to Aq without making any more progress (t^). The term in the denominator is for normalization 
and represents the probability of returning to A. The expression for Tf is much simpler: Tf — X]j=o ^l^^ ■ These 
values oiTint, Text, and Tf can be combined with P(Am|Ao) in (fll) to obtain Ttot- By measuring Text in the second step 
of FFS, the value of Ttot obtained will be comparatively more accurate and have less variance. The only extra work 
done in FFST over FFS is keeping track of the simulation time during the second step, which makes the advantages 
gained by FFST practically free. 

B. Forward Flux with Quasi-Markov Dynamics 

In order to gain more insight into the connection between the FFS and the barrier method, we formulate FFS as 
quasi-Markov dynamics. That is, we iteratively calculate the time it takes to travel between three barriers. The first 
step is similar to FFS: run a single long simulation and calculate the time it takes to reach the first barrier (Ai) from 
the surface of A (Aq), and where along Aq the sample crosses going out. N samples are started along Ai, where the 
initial path crossed in the first step, and are run until they reach A2 or Aq. The average times to go from Aq to Ai, Ai 
to A2, and from Ai to Aq are tj, if, and i^, respectively. In general, the time takes to reach A^ from Ai is given by tj. 
The probability of reaching A2 from Ai without first going back to Aq is P(A2|Ai). These times and probabilities can 
be used to make a random walk with three states, Aq, Ai, and A2. The transition time from Aq to A2 (tp) is given by 
the weighted times of all possible paths from Ao to A2 . These paths can be organized by the number of times they 
return to Aq. Writing out the first few terms in this series exposes the general form, 

tl = [{tl + tl)P{X2\X^)] + [{tl +tl + (i? + 4))P(A2|Ai)(l - P(A2|Ai))] 
+ [{tl + tl + 2(t? + ii))P(A2|Ai)(l - P(A2|Ai))2] + . . . 

tl + k{tl + ti))P(A2|Ai)(l ~ P(A2|Ai))'=. 



= E(*o 



fc=0 



The above expression is the sum over all possible ways to reach A2 from Aq, sorted by the number of times the 
simulation returned to Aq. The probability of returning to Aq k times before reaching A2 is given by the product of 
the probability of not reaching A2 k times ((1 — P(A2|Ai))'^) and then reaching A2 on the {k + 1)*^ try (P(A2|Ai)). 
The time this takes is given by the sum of the time it takes to go from Aq to Ai to Aq k times {k{ti + 1\)) plus the 
time it takes to make it from Aq directly to A2 {t\ +t\). Note that in practice the sum converges quickly because of 
the factor, (1 - P(A2|Ai))'=. 

The next step is to repeat the same process using Aq, A2, and A3 as the three barriers. Brute- force dynamics are 
used to measure t^, ^2' and P(A3|A2). Then, the three-barrier calculation from Q is used with 0, 2, and 3, in place 
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of 0, 1, and 2. The result from the calculation is an estimated value of tg. In general to calculate t^Q^^ for the i 
we use: 



th 



step 



A+i 



oo 

E 

k=Q 



[fo + tf^ + k{t\ + ij,)] P(A.+i|AO(l - P(A.+i|A,))' 



(5) 



This process of performing short brute- force simulations, followed by solving ([5|, is repeated for every barrier until A 



M 



is reached. At the end we have t^ — Ttot- This method could be useful for practical simulations. Here we introduce 
it as a pedagogical device to show that by making small changes to the FFS algorithm the barrier method can be 
effectively obtained. 



Forward Flux Barriers 



We can look at quasi-Markov dynamics in another way. We start by measuring the average time it takes to reach 
Ai starting at Aq (ij), while also keeping track of where along Ai the path crosses going out. The paths are continued 
from Ai until they reach A2 or Aq. If a sample reaches Aq then it is restarted at Ai at one of the locations where paths 
ended in the first step, and t\ is added to the time. This process is continued until all samples reach A2. We now 
have an estimate of t^ and the locations along A2 where the sample paths ended. From these locations the paths are 
continued until they reach A3 or Aq. If they reach Aq they are restarted at a location where a previous path stopped 

This step is finished once all sample paths reach A3. The general step is to start the 

with in added to the 



A2 and t\ is added to the time 



This process is repeated until Am is 



paths on Afc and run them until they reach Afc+i or Aq. If a path reaches Ao, it is restarted at A^ 
time. The step is complete when all paths reach Afc+i and time time gives ig^^. 
reached. The result is a value for t^ which is the estimate for the transition time. 

This construction shows the relationship between FFS and the barrier method |36j . This version of FFS measures 
the average time it takes to reach each barrier during the algorithm, as in the barrier method. Also both methods 
avoid characterizing the flux rate from the surface and the transition probability. The main difference between this 
algorithm and the barrier method is that the barrier method need the simulation to go all the way back to Aq before 
jumping back to the current barrier; only the previous barrier need be reached. This is the source of the performance 
gains of the barrier method. However, this method is currently only tractable for low-dimensional systems as it 
requires a reasonable sample of previous barriers. 



VII. APPENDIX B 



TESTING FFS AND FFST ON AN EXACTLY SOLVABLE PROBLEM 



In this Appendix, we use a simple one-dimensional system to show that FFST can give accurate results for transition 
times when FFS fails. We also briefly discuss the comparative efficiency of the different algorithms on the ZGB model. 

Consider a discrete hopping process on a line of length L. The probability of jumping from the i*^ to the [i — 1)'*' 
site is "Pi. The time of a jump is unity. In the cases where the first site (i — 0) is adsorbing and the last site is 
refiecting (j)h-\ = 1), the system can be solved exactly |40j . The solution can be written in terms of hopping rates 
instead of hopping probability. Ai and /ij are the rates of hopping from i to (z -I- 1) and (i — 1) respectively. In terms 
of "Pi these are: \Xi = pi and Ai = I — Pt- The average time to reach the adsorbing site i = from site i — n is [40] : 



E 

m—l 

n 

E 



— J- Y\ — \^ J_ TT ^ 

m-l L-1 1 J-1 / 1 

P^ f=i l-PS=»+l^■'fe=A^^■ 



(6) 



This equation has been used to find the extinction time of a disease within a population in a simple model from 
epidemiology ^40]. Equation ^ is general and can be used to construct 'energy landscapes.' We create a landscape 
with non-uniform hopping probabilities such that there are three metastable states, regions A, B, and C, all with 
roughly equal stability, as shown in Fig. [TO^ . We measure the time it takes to reach the absorbing state starting 
near the reflecting boundary. This requires escape from the first metastable region A, then the second metastable 
region B, to finally reach the absorbing state near the center of the last metastable region C We measured the time 
using FFS and FFST for various sizes of systems (well depths) as shown in Fig. 10 3. We found that FFS significantly 
underestimates the transition time by as much as 50%. There is also a significant increase in the variance of the 
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result, roughly an order of magnitude for this model system. This is caused by the flux being strongly influenced by 
the rare occurrence of a trajectory that makes it to region B, spends a long time there, and then returns to region A. 
Even in the absence of long-lived metastable states, FFST can produce transition times with less variance than 
FFS because it samples the external return time signiflcantly better. We found this to be the case in the ZGB model. 
Using L — 32, eleven evenly spaced static barriers starting at C = 0.06, N = 10'^, and pco = 0.5268, we found a 
variance of 5.0% for FFST and 5.8% for FFS; an improvement of about 16%. To equal the variance of the FFST 
result, FFS which would require roughly 35% more trials which would translate into a 35% longer run time. In the 
case of the ZGB model, this effect can be mitigated by choosing Aq to be far enough from the metastable region that 
the internal return time Tint is much larger than the external return time Te^t, bounding the effect on the increase 
in variance of Text- In general, moving Aq is not always useful, because metastable states can make Text arbitrarily 
large. 
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